Distributional exact diagonalization formalism for quantum impurity models 
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We develop a method for calculating the self-energy of a quantum impurity coupled to a continuous 
bath by stochastically generating a distribution of finite Anderson models that are solved by exact 
diagonalization, using the noninteracting local spectral function as a probability distribution for the 
sampling. The method enables calculation of the full analytic self-energy and single-particle Green's 
function in the complex frequency plane, without analytic continuation, and can be used for both 
finite and zero temperature at arbitrary fillings. Results are in good agreement with imaginary 
frequency data from continuous-time quantum Monte Carlo calculations for the single impurity 
Anderson model and the two-orbital Hubbard model within dynamical mean field theory (DMFT) 
as well as real frequency data for self energy of the single band Hubbard model within DMFT using 
numerical renormalization group. The method should be applicable to a wide range of quantum 
impurity models and particularly useful when high-precision real frequency results are sought. 



PACS numbers: 74.25. Ha,74.25.Jb,74.72.-h,79.60.-i 

The topic of quantum impurity models is a corner- 
stone in the study of strongly correlated electrons. The 
standard model is the Anderson impurity model which 
was first developed to describe dilute magnetic ions in 
metals. It contains the most classic manifestation of 
strong correlations, the Kondo effect. A modern appli- 
cation of quantum impurity models has emerged in the 
study of strongly interacting lattice models using dy- 
namical mean field theory (DMFT) EE where the peri- 
odic system is mapped to an quantum impurity model 
coupled to a non-interacting bath, by assuming that cor- 
relations described by the self-energy, are fully dynamic 
but local. Although this assumption is exact only in the 
limit of infinite dimensions, DMFT has proven of great 
value in the study of finite dimensional correlated sys- 
tems. Despite of providing a drastic simplification com- 
pared to the full momentum dependent lattice problem, 
solving the quantum impurity problem is still a major 
challenge. 

With the advent of high resolution angle resolved pho- 
toemission spectroscopy as one of the main probes of 
strongly correlated systems there is an increasing de- 
mand for calculations of the single particle spectral func- 
tion, using for example DMFT. In this context the most 
versatile and widely used impurity solver is the contin- 
uous time quantum Monte Carlo (CT-QMC) methodP 
The sampling of the thermal Green's function in imag- 
inary time or frequency is very efficient but, statistical 
errors grow with decreased temperature. Calculating 
the spectral function from the CT-QMC results requires 
a numerical analytic continuation to the real frequency 
axis. Here the statistical errors prevents the use of Pade 
based meth ods ^ while the maximum entropy method 
(MAXENTpE can be used to find a smoothed out spec- 
tral function. To be able to capture the detailed spec- 
tral function it is clearly preferable to work directly with 
real frequencies or time. The numerical renormalisation 
group formalism (NRG does exactly this and has been 
used to extract fine details such as small kinks in the 
quasiparticle dispersion caused by correlation effectsP^Sl 
Nevertheless, the method is very demanding for models 
with more than one impurity site. 

Another very versatile tool for DMFT calculations is 



the so called exact diagonalization (ED) formalism^ in 
which the quantum impurity problem is represented by 
an exactly solvable finite Anderson impurity model with 
a small set of non-interacting bath levels. The fit of the 
few bath levels to the continuum non-interacting impu- 
rity Green's function can be surprisingly accurate on the 
imaginary frequency axis and it is possible to calculate 
the corresponding self-energy to very good accuracy.^ 
However, if the impurity self-energy is evaluated along 
the real frequency axis the finite-size effects of the few- 
level Anderson model become apparent, giving a dis- 
crete set of poles (see eg. Fig. 18 in Ref. [T]), while the 
self-energy of the impurity coupled to a continuous bath 
should be continuous. 

In this paper we present a formalism for calculating 
the full analytic (real and imaginary frequency) self en- 
ergy of a quantum impurity model by using a stochastic 
sampling of the non-interacting impurity spectral func- 
tion instead of the Matsubara fit of the ED formalism. 
Each sample is a small n- level Anderson impurity model 
for which the self energy can be calculated by exact di- 
agonalization. By sampling over the impurity Green's 
function the full phase space can be explored and a, for 
practical purposes, continuous self energy is formed as 
a sample average. 

In terms of applications we believe that the method is 
as versatile as the ED formalism but with greater phys- 
ical relevance. Comparisons with CT-QMC shows very 
good agreement on the Matsubara frequencies as well 
as with NRG real frequency results for the self energy. 
Compared to CT-QMC a major strength of the method 
is that real frequency results are calculated without any 
analytic continuation and that the method is most ef- 
ficient for zero temperature. Because a large stochas- 
tic sampling is required for accurate results the calcu- 
lations can be numerically demanding. However, the 
samples can be generated and addressed independently 
making it ideal for parallel computing. In this paper 
which introduces the method we will focus on the single- 
impurity Anderson model with a semicircular bare den- 
sity of states, and the corresponding single-orbital Hub- 
bard model when considering the application to DMFT, 
but also show a DMFT calculation for a two band Hub- 
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FIG. 1: Basic approximation of the formalism corre- 
sponding to grouping of poles (Eq. [T| and ignoring cross 
correlations between groups, illustrated here for a 2nd order 
diagram. 

bard model. 

For the single-impurity Anderson model the basic 
object of study is the imaginary time action S = 
-Sdrdr'^cUr^ir-r') + »5(t t')]c ct (t') + 
U J drcl(r)cf (r)c4.(r)c-|-(r), with spine =t, 4- The non- 
interacting impurity Green's function Go describes the 
correlations induced by the coupling to the surround- 
ing non-interacting bath. In complex frequency space, 
Gq(z) is an analytic function with poles (branch cut) 
on the real frequency axis. The task is to calculate 
the Green's function G a {r - r') = -(Tc ct (t)4(t'))s. 
Subsequently we will assume no magnetic order and 
drop the spin index and instead of the Green's function 
we can consider the self energy £ given by G" 1 (z) = 
Gq 1 {z) + /j, — For the quantum impurity prob- 

lem, Go corresponds to the bare (non-interacting) den- 
sity of states which we take to be semicircular Po(uj) = 
(2/7r)vI — w 2 and let the half-bandwidth be our unit of 
energy. 

Consider a representation of Go in terms of a large 
number of poles on the real axis Gq(z) = X)f=i 
with ^ 5j = 1, which for the purpose of doing numerical 
calculations can be very good for a large number (10 4 or 
more) of poles. Assume that the poles and residues are 
such that they can be grouped in N smaller groups of 
size n (the total number of poles is thus M = Nn) such 
that the total residue in each group is 1/N. (There are 
many inequivalent ways of grouping the poles, we make 
an unbiased choice, grouping the poles randomly.) We 
rewrite 
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where the residues are renormalized by a factor N such 
that the Green's functions Gq(z) are properly normal- 
ized for an n-level system, ^ ■ a J = 1. 

The self energy is given by all one particle irreducible 
(1PI) diagrams in terms of the four point vertex U and 
the two point vertex —fi connected by Go-^ For every 
diagram, we make the following approximation exem- 
plified in Fig. [I] by a 2nd order diagram (omitting fx 
insertions) . 

The approximation is thus to replace cross correla- 
tions between different Green's functions by internal 
correlations, giving £ « A J2 U ^ where is the self 
energy related to Gq. 

Considering also the chemical potential \x on the in- 
teracting site, all 1PI diagrams include the diagram 
with the vertex —\x by itself as well as all insertions 
of —fx into the 1PI diagrams constructed from the 
vertex U. Within the same approximation, /i acts 
as a chemical potential on each subsystem such that 
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find that the approximation corresponds to the expres- 
sion 
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which is the basis of the present formalism. 

Importantly, E" contains all 1PI diagrams of Gq, it is 
the exact self energy of the quantum impurity action S, 
with Go replaced by Gq, a problem that can be mapped 
to an Anderson impurity model with a single interacting 
site coupled to n — 1 bath levels. The Anderson model 
is formulated in terms of a Hamiltonian which can be 
diagonalized numerically for small n and the self energy 
calculated as 



^{z)-fi={G^z))- 1 -{G v {z)r 



(3) 



Note that it is not an option to work with the sample 
averaged Green's functions instead of the self energy; 
jj Y^, v G v and J2 U Gq is not a proper pair of interact- 
ing and non-interacting Green's functions, they do not 
obey Go(z) = =>■ G(z) = 0, as required by the Dyson 
equation. 

Now, Gq = X)j=i z-b v * s ma PP e d to the Green's func- 
tion Gl = l/(z - e v Q - YTjZl ^r) of the Anderson 
model H = eg £ CT c\c a + E"J=iV/(4ci<r + h.c) + 
ejCj a Cj a ] by solving for the parameters and Vi ac- 
cording to et : Gg(w = e\) = 0, ^ = -l/(^) 2 , 
and t v o = -T,j a j b j' The ful1 Hamiltonian is H = 

Ho ~ A 1 X)o- c o- c ct + ^ c t c l c ^ c T ano - the corresponding 
Green's function G v (z) given by the Lehmann repre- 
sentation by summing over the complete set of eigen- 
statesP 

The n-level systems is derived from a representation 
of the full Go in terms of a large number of poles 
grouped into sets of poles with normalized residues, 
Eq. [I] An exact representation of Gq(z) is a distri- 
bution of poles given by the spectral function Ao{uj) — 
— ^Im[Go(w + i0 + )] and in practice we will use this as 
a probability distribution for generating sequences of 
n random pole locations with equal residues. (Other 
sampling procedures are conceivable.) A large ensem- 
ble of such groups will give a good representation of 
Go- However, one show that the approximation (illus- 
trated in Fig. [T]) overestimates the contribution of low 
energy spectral weight in Go- We have found that this 
can be compensated by imposing the constraint that 
the ground state particle number of the interacting and 
non-interacting n-level systems coincide, discarding con- 
figurations that fail to satisfy this criterion. (This also 
means that the sampling will not be a completely faith- 
ful, as seen in Fig. |4j) 

The suggested operational procedure for a particle- 
hole symmetric system is: 

• Use Aq{u) as a probability distribution for generating 
N sets of n-poles, with normalized residues, giving the 
Green's functions Gq, v = 1, N. 

• For each v check that the ground state of the inter- 
acting Hamiltonian has the same particle number as the 
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FIG. 2: (Color online) DMFT real frequency self-energy 
E(w) (solid line) for U=l, T=0, at and away from half- 
filling, fi = 0.5 and n — 0.75 respectively, corresponding 
NRG results (dashed line) from Crete et al, Ref.QJJ Model 
size n = 6 and 2 • 10 4 samples with acceptance ratios 89% 
and 81%. 

non-interacting, else discard the configuration v. Cal- 
culate the self energy E" by mapping the problem to a 
Hamiltonian and using Eq. [3j 
• Add up the self energy E = ^ Y, v s " 

When applied to DMFT Go is calculated from the 
DMFT equations for the local Green's function Gl 

Gl{z) = /f ^_J?ff s(z) and G- Q \z) = G7» + 
Y<(z) — fi which follows from integrating out all degrees 
of freedom except those of one single (impurity) site. 
Given Go the impurity problem is solved for the self- 
energy E which gives a closed set of equations that are 
solved self-consistently. 

Away from half filling the standard formalism is 
poorly suited for the sampling, as the particle number 
of Go and Gl may be very different. We will use a for- 
mally equivalent expression, defining E(z) — E(z) — So 
where Eq is a real constant. In terms of this we write 



G L (z)= (Go(z))- 1 - E(z) 



(4) 



where G 1 (z) = G 1 (z) + /i — E and choose E such 
that J du> lm[Go(uj)]nF(oj) = J dujlm[Gi,(Lj)]nF(^)-> the 
two Green's functions give the same occupation. Thus 
starting each step of the DMFT iteration we would solve 
for Eo by fitting the particle numbers, use Eo as the 
effective chemical potential for the n-level systems and 
the spectral function of Go as probability distribution. 
(For half filling E = /i = u/2 and this step is trivial.) 
For the quantum impurity model we work with G = 
l/((Go,fc are ) _1 + A 1 — instead of Gl and choose E 
such that G and Go = l/(Gg 1 +/x— E ) have the same 
particle number. Here we need to converge E , because 
E enters into G and the particle numbers of G and G 
should be same for the best sampling. In what follows, 
we present several different calculations and compare to 
NGR and CT-QMC. 

In Figure [2] we show DMFT calculations for U = 1, 
T = 0, with [i = 0.5 (half filling) and fi = 0.75 to 
compare with published real frequency NRG- results 
for the self energy. For the system away from half filling 
we use a sampling where the particle number of the 




FIG. 3: (Color online) DMFT calculation for U = 2.8, 
T — 0, (upper left panel) real part of S(oj) — U/2 (solid line) 
with kinks indicated fitted dotted lines, (upper right panel) 
full view of the Self-energy, (lower panel) —lm[G(z)} in the 
upper complex plane (grid-surface) and on the real axis 
(solid line), coexisting insulating solution (thin solid line). 
Model size n = 6 and 3 • 10 4 samples with 49% acceptance 
ratio. 

non-interacting samples are free to vary between and 
2n (in units of 2) guided by the probability distribution. 
The results are in strikingly good agreement considering 
that there are no free parameters in the formalism. 

In Figure [3] we show results for U — 2.8, close to 
the T = Mott transition at U w 3, 8 both the metal- 
lic and insulating solutions. Clearly, the method is 
sophisticated enough to capture fine structure in the 
self energ y, co rresponding to kinks in the quasiparticle 
dispersion!21i2l 

Figure [4] is a finite temperature calculation at j3 = 
30 (T = 1/30) for the Anderson impurity model at 
U = 3. We find good agreement with CT-QMC cal- 
culations for the impurity Green's function G{z) — 
1/(^0 lare( z ) +M ~ ^( z )) on the Matsubara frequencies 
z = iuj n = i^(2n + 1), with maximum error for n = 6 
is 4%. Interestingly, the real frequency density of states 
is quite different from that deduced from MAXENTpl 
The central (Kondo) peak is similar but the Hubbard 
bands are much more distinct. 

To show the generality of the method we also present 
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FIG. 4: (Color online) Anderson impurity at half-filling, 
[7 = 3 and j3 — 30. Interacting Green's function G(z) on the 
imaginary axis (solid line) compared to CT-QMC results 
(crosses) using the Triqs-codep^Sl bath Green's function 
Go (dashed line). (inset) p(t-u) = — Im[l/7rG(cj + iS)] 
(solid line), bath DOS po(^) (dotted line), and sampled po 
discarding configurations (dashed line). Model size n — 6 
and 10 4 samples with 65% acceptance ratio. 
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FIG. 5: (Color online) Two band Hubbard model at 
T = with U = 1, J = 0.2, and half-band widths Wi = 1 
and Wi = 0.1 in the orbital selective Mott phase. With 
density-density (dashed lines) and full Hund's (solid 
lines) interaction, the former is a non-Fermi liquid with 
Im[Ei(0)] 7^ 0. (Top) Im[G(z)] along the imaginary axis 
compared to CT-QMCP™ (j3 = 120) for the wide (left) 
and narrow (right) band. (Bottom) Real-frequency spectral 
functions (narrow band scaled by 1/10), non-interacting 
case (dash-dotted lines) density-density (dashed lines), and 
Hund's (solid lines) interaction. Model size n — 6 and « 10 4 
accepted samples. 

a calculation for the two band Hubbard model both with 
density-density and full Hund's rule coupling, for model 
details see Refs. 1201211 The spectral functions for the 
bands in the orbitally selective Mott phase are shown 
in Fig. [5] and the imaginary frequency Green's functions 



are compared to CT-QMC results. The sampled sys- 
tems have a total size of n — 6 with three (one) bath 
sites for the metallic (insulating) band. For the sys- 
tem with only density-density interactions, which is a 
non-Fermi liquid, we have used the fact that the spin is 
conserved in each orbital, to calculate two independent 
self-energies for the wide band, depending on the spin 
in the narrow band, see Biermann et al. Rcf. 1211 

In summary, we have presented a formalism for calcu- 
lating the full analytic self energy of quantum impurity 
models by using a representative distribution of exactly 
solvable Anderson impurity models. The method is sim- 
ple to implement and the initial studies shows that the 
method can give very good results. The calculations 
in this paper were done a single desktop computer over 
time periods of 10-40 hours, but the formalism is well 
suited for parallel computing which will be the key to 
considering larger n models. A natural extension is to 
apply the formalism to more general models, including 
multi orbital and cluster generalizations of DMFT and 
to calculate other dynamical correlations. 
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